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I.  Introduction 

The  main  problem  of  sound  ranging  is  to  determine  the  location  of  a 
sound  source  by  examining  the  signals  received  at  an  array  of  microphones. 
For  example  consider  the  two  dimensional  problem  with  a sound  source  on  the 
same  plane  as  a linear  array  of  six  microphones  as  in  Figure  1. 
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Figure  1.  Linear  Array  Sound  Ranging 
The  sound  wave  travels  outward  from  the  source  having  roughly  a spherical 
wavefront.  The  time  the  wave  takes  to  reach  each  microphone  depends  upon 
the  temperature,  atmosphere,  wind  profile  and  distance  between  the  source 
and  the  microphone  as  well  as  the  terrain  if  the  array  were  positioned  on 
the  ground.  In  Figure  2 an  example  of  the  signals  from  the  outputs  of  the 
configuration  in  Figure  1 for  a transient  source  are  given.  These  signals 
are  typical  of  those  received  from  a Howitzer  or  explosion  excluding  any 
ballistic  waves. 
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Figure  2.  Microphone  Signals 
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Most  existing  sound  ranging  solutions,  for  example  Field  Manual  [1], 
Miller-Engebos  [ 2 ],  and  Bangs-Schultheiss  [ 3 ],  use  the  time  differences 
between  the  sound  bursts  and  the  existing  temperature  and  wind  information 
at  the  microphone  sites  as  critical  information.  It  is  thought  that  better 
estimates  of  the  time  differences  between  bursts  will  provide  better  esti- 
mates of  the  source's  location.  Many  existing  techniques  determine  the 
time  difference  by  first  selecting  individual  "break  points"  on  each  of 
the  microphone  signals  and  then  subtracting  to  obtain  the  time  differences. 
Some  of  the  methods  used  to  determine  these  break  points  are; first  inflec- 
tion from  the  noise  baseline,  first  maximum  after  inflection, or  first  zero 
cross  over  after  inflection.  In  contrast  to  these  methods  the  general 
approach  presented  in  this  paper  is  one  of  estimation  of  time  differences 
between  signals  by  using  the  normalized  cross  correlation  functions  between 
filtered  versions  of  the  signals.  In  this  way  the  entire  structure  of  the 
signals  is  used  rather  than  a sometimes  vague  "break-point"  and  error  caused 
by  noise  can  be  minimized.  After  these  time  differences  have  been  obtained 
the  location  of  the  sound  source  can  be  determined  by  any  of  the  existing 
algorithms,  some  of  which  are  compared  in  Engebos  and  Miller  [4]. 

II.  Mathematical  Model 

In  order  to  estimate  the  time  differences  a mathematical  frame  work 
must  be  established. 

The  general  model  to  be  used  in  this  paper  is  shown  in  Figure  3.  The 
sound  source  is  assumed  to  generate  a compression  wave  that  is  de- 
fined as  p(t).  As  the  wave  travels  through  the  atmosphere  or  along  or 
through  the  ground  it  is  changed  as  if  it  has  passed  through  a filter 


of  some  kind.  Although  in  general  this  filtering  action  is  time  varying 
and  distributive  we  will  assume  a linear  time  invariant  filter  with  fre- 


j 


quency  response  A^co)  followed  by  a time  delay  t^.  The  subscript  indicates 
filtering  of  the  wave  before  it  reaches  the  iC^  microphone.  Each  micro- 
phone itself  has  a frequency  response  M^Cui).  The  microphones  when  positioned 
sufficiently  far  apart  each  receive  independent  background  noise,  wind  noise, 
and  other  extraneous  noise. 

It  is  assumed  that  these  noise  sources  can  be  collectively  modeled  as 

wide  sense  stationary  independent  random  processes  r|^(t)  with  power  spectral 

densities  <J>..(oj).  The  basic  problem  then  is  to  estimate  the  times  t.  after 
ii  1 

receiving  the  signals  x^Ct):  i=l,2,...,6  over  the  time  interval  from  tr  to  t^. 
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The  general  model  can  be  simplified  to  that  shown  in  Figure  4 if  we 
assume : 

(1)  A^io)  = a^Atco),  that  is  we  have  identical  filtering  in  all  paths 


but  a different  attenuation  due  to  the  distance  traveled  from  the 


source  to  each  microphone. 

(2)  M^(oo)  = M(w),  that  is  all  microphones  have  the  same  frequency  response. 


(3)  n.(t)  A n. (t)  * m(t)  represents  the  background  noise  n.(t)  reflected 

1 - l i 


through  the  microphone  where  m(t)  is  the  impulsive  response  of  the  micro- 
phones and  * means  convolution  in  the  time  domain. 

(4)  q(t)  A p(t)  * a(t)  * m(t)  represents  the  undelayed  signal  version  in- 
cluding the  filtering  action  of  the  atmosphere  and  microphone.  The  a(t) 
is  the  impulsive  response  of  this  filtering  action,  i.e.,  a(t)  =?  1 [ A (uo)  ] 


q(t) 


x (t)«a  q(t-t  )+n (t) 


h (c) 

~*\+)r x2(t)=a2q(t-t2)+n2(t) 

x3(t)=a3q(t-t3)+n3(t) 


x^  (t)=£*^q  (t-t^)-Hi^  (t) 


x5(t)=a5q(t-t5)+n5(t) 


Jn  (t) 

-© * x6(t)=abq(t-t6)+n6(t) 


Figure  4.  Simplified  Model 

For  any  two  of  the  signals  x^t),  *j(t)  where  i t j , and  the  following 


definitions  for  s(t),  B„  and  T„ 


s(t)  = Ctiq  (t-t ±)  , 


6 A ^i_  A 

6ij  ' a.  ’ Tij  ‘ Vj 


(1) 


we  have 


x3(t)  = s ( t ) + nt(t) 


(2) 


x j ( t ) = 8ij  s ( t— T 3 j ) + nj(t) 
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The  problem  now  becomes  the  estimation  of  . from  the  received  sig- 
nals x^(t)  and  Xj  (t) . This  simplified  model  will  be  used  as  the  basic 


model  for  the  paper.  Hannan  and  Thompson  [ 5 ] have  developed  a maximum  likeli- 


hood estimator  for  T and  Knapp  and  Carter  [6]  present  several  other  estima- 

ij 

tor  forms  for  various  performance  characteristics.  The  basic  structure  of 
these  estimators  has  been  shown  to  be  a combination  of  prefiltering  and  cross- 

A 

correlation  as  given  in  Figure  3.  The  time  argument  at  which  the  cor- 

A 

relator  reaches  a maximum  y . is  the  delay  constant.  The  t and  y for  all  i 

ij  i j i j 

and  j thus  determined  can  then  be  used  in  any  of  the  existing  algorithms  to 
estimate  the  location  of  the  sound  source. 

' 1 I 1 
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Figure  5.  Structure  of  Estimator  for  a Pair  of  Signals 
The  prefilters  are  determined  by  the  type  of  estimator  used  and  assump- 
tions on  the  noise  source  and  sound  source;  examples  are  given  by  Carter  and 
Knapp  [ 6 ] . In  the  usual  case  the  actual  frequency  content  of  the  signals 
may  be  unknown  as  well  as  the  constants  and  the  noise  power  spectral 
densities.  Forms  of  filters  H^(jo))  and  H2(joj)  might  then  be  formed  using  on 
line  estimates  of  these  parameters.  Estimates  of  the  noise  spectral  densities 
of  each  signal  can  be  easily  obtained  prior  to  the  signal  epoch  if  enough 
lead  time  is  given  by  using  an  FFT  algorithm.  This  ^oise  energy  level  together 
with  the  signal  pulse  noise  energy  could  be  used  to  estimate  a and  a filter 
H^(ju)  and  H2 ( j w)  formed  from  these  estimates. 
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III.  Proposed  Method 


From  the  six  microphone  signals  the  location  of  a sound  source  is 
desired.  The  method  proposed  for  determining  the  sound  source  location  is 
a digital  processor  composed  of  two  basic  parts.  The  purpose  of  the  first 


is  to  estimate  the  appropriate  time  differences  between  the  arrival  of  sound 
source  signals,  while  the  purpose  of  the  second  is  to  use  these  estimates 
to  obtain  an  estimate  of  the  source  location.  The  first  section  is  the  main 
subject  of  this  report  while  details  of  the  implementation  of  the  second  on 
a desk  computer  (HP  9825)  are  given  by  Miller  and  Engebos  [ 7 ]. 
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Figure  6.  Digital  Processor  Structure 

The  time  differences  estimator  used  is  a combination  of  optimal  estima- 
tor structures  (Figure  5)  for  all  possible  pairs  of  signals  followed  by  selec- 
tion of  one  of  the  signals  as  a reference  (call  it  i^  e[l,2 , . . . ,6] . The  over- 
all structure  of  the  time  differences  estimator  is  given  in  Figure  7 and  the 


details  are  presented  in  the  following  section. 
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Figure  7.  Time  Difference  Estimator 
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IV.  Time  Difference  Estimator 


The  time  difference  estimator  can  be  broken  into  three  subsections  :(A) 
Estimation  of  rough  time  differences  and  measure  of  correlation  for  all  pairs 
of  signals  (B)  Elimination  of  unreliable  time  difference  information  (C)  Re- 
adjustment of  time  difference  estimates  for  overall  best  fit.  The  details 
of  each  of  these  subsections  follow. 

A.  Rough  Estimation  of  Time  Differences 

Each  of  the  received  signals  x^(t)  has  been  received  over  the  time  in- 
terval t to  tj.  In  many  cases  the  signal  duration  that  is  being  searched 
for  is  a small  percentage  of  this  total  interval.  The  structure  shown  in 
Figure  5 could  be  used  throughout  the  interval  from  to  t^.  but  the  actual 

correlation  done, either  direct  or  with  the  Fast  Fourier  Transform, may  take 
an  exorbitant  amount  of  time.  Since  the  signal  to  be  located  is  a transient 
type  it  is  more  efficient  to  determine  rough  starting  times  for  the  source 
signals  by  using  an  energy  detector  and  then  use  the  cross  correlation  as 
a vernier.  The  overall  block  diagram  for  the  estimation  of  rough  arrival 
times  t^  is  given  in  Figure  8. 
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Figure  8.  Rough  Starting  Time  Estimator 
The  signal  Xj(t)  or  the  output  of  Filter  F^  from  x^(t)  Is  first 
squared  then  integrated  over  a sliding  window  of  length  Ty  beginning 
at  t and  ending  at  t+Ty  where  t Is  varied  from  tr  to  t^-Ty.  The  Ty 
Is  selected  to  be  approximately  the  same  length  as  the  s(t)  expected. 
The  energy  detector  can  be  operating  on  the  data  in  real  time  and 

A 

the  time  at  which  a maximum  Is  reached  gives  a rough  estimate  tj 
of  the  relative  time  of  arrival  of  the  signal  we  are  looking  for  at 
each  of  the  microphones. 

Operating  on  all  of  the  six  signals  In  this  fashion  produces  a set 

A A A 

of  rough  arrival  times  t.,  t_,...,t,.  In  this  wav  we  can  narrow  the 

t Z o 

cross  correlation  operation  window  and  reduce  the  calculation  time 
for  estimating  the  rough  time  differences.  The  normalized  cross  cor- 
relation for  a pair  of  signals  x ^ ( t ) and  x^ft)  is  then  given  by 

V  l * * *  5V‘ 


VT)  ' 


Xj(t)  x | (t  + i )dt/  /k j K 


(3) 


ll  ‘w7^  where  ■ 


V 3Tw/2 
x.‘ (t )dt 


*i-  V2 


The  total  time  of  the  Integration  Is  selected  as  3/2  Ty.  This  makes 

sure  the  entire  signal  will  be  within  the  cross  correlation  window 

provided  the  time  of  arrival  Is  roughly  guessed  to  within  one  quarter 

of  its  duration  in  either  direction. 
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If  t <t  then  T Is  varied  from  t.-t.  -T.,/2  to  t.-t,  + T /2  while  if 
ji  ijW  ijW 

t >t  , T is  varied  from  t.-t  - Tu/2  to  t,-t  +Tu/2.  The  value  T fo,r 
ji  jiwjiw  ij 

which  Y^j(t)  reaches  its  maximum  y^j  is  a rough  estimate  of  the  time 
difference  between  the  signals  x^(t)  and  x^ (t)  , while  the  Yjj  gives 
a measure  of  the  correlation  between  the  signals.  The  Y^j(t)  can  be 
computed  directly  by  time  integration  for  various  values  of  t or  by 


the  use  of  the  FFT. 

B.  Elimination  of  Unrealiable  Time  Difference  Information 

In  some  cases  there  may  be  very  weak  correlation  between  the  sig- 
nals x (t)  and  Xj  (t) . This  weak  correlation  could  be  caused  by  sig- 
nal s(t)  being  very  small  compared  to  the  additive  noise  or  by  the 
fact  that  s(t)  has  been  severely  distorted  in  transmission  shedding 
doubt  on  the  simplified  model  of  Figure  A.  In  either 

case  the  rough  time  difference  estimate  is  inaccurate,  sometimes  to  the 
point  of  being  totally  unreliable.  If  these  unreliable  time  differences 
were  used  in  the  final  refinement  or  readjustment  explained  in  part  C, 
severe  errors  in  the  refinement  may  result  and  consequently  be  the 
cause  of  large  errors  in  the  desired  range  and  bearing  estimates. 

If  y^j  is  small  compared  to  one  it  means  that  the  signals  x^Ct) 

/v 

and  Xj (t)  are  weakly  linearly  related  and  the  time  difference  is 
unreliable  while  if  Yjj  is  approximately  one,  the  signals  are  strongly 
linearly  related  and  the  time  difference  estimate  is  reliable.  Define 


a threshold  T^  in  such  a way  that  Y^j  < T^  would  imply  the  time  dif- 
ference estimate  is  unreliable  and  Y^^  Tq  implies  a reliable  estimate. 

In  this  way  a signal  x^(t),  1*1, 2,.., 6 ^or  which  Y^  < for  all  j=l,2,...6 

A 

would  be  eliminated  from  consideration.  If  Y^j < Tq  for  just  some  j's, 
but  not  all  of  them,  those  T^'s  should  be  eliminated  from  consideration 
but  not  the  entire  signal. 
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B.  • 


C.  Readjustment  of  Time  Difference  Estimates 

After  the  unrealiable  time  differences  have  been  eliminated  a best 
overall  fit  may  be  obtained  by  a least  square  procedure.  The  following 
presentation  is  for  when  none  are  eliminated.  If  elimination  occurs,  a 


renumbering  must  take  place.  By  selecting  one  of  the  signals,  say  x.,(t) 


as  a reference  (if  another  is  chosen  simply  renumber)  the  time  difference 


of  signal  x^(t)  relative  to  x^(t)  for  i = 2, 3,..., 6 can  be  calculated  as 


(see  appendix  for  derivation) 


6 (T  + t12  - x23 


T24  ~ T25  ~ T26] 


6 [T  + T13  + X23 


34 


35 


T36] 


6 [T  + t14  + t24 


4-  T - T 

34  45 


T46] 


(5) 


6 + T15  + T25 


+ T35  + T45 


T56] 


6 [T  + T16  + T26 


+ T + T + T 1 

36  46  56J 


where 


'■'l!*  Tn  + Tu  + T15  + t16 


(6) 


The  reference  signal  x (t)  chosen  should  be  the  one  for  which  the 

R 


sum  of  the  y . for  all  j^R  is  a maximum.  If  some  of  the  data  is  un- 

Rj 


reliable  we  may  proceed  as  follows.  Suppose  is  less  than  T^  thus 


giving  rise  to  an  unreliable  If  we  8°  ahead  and  use  this  un- 


reliable estimation  of  the  above  formulas  we  adversely  affect  our  esti- 
mates. Nonetheless  we  must  still  substitute  something  into  equations 


for  T_ . . If  we  can  find  a k such  that  yT1  and  y.,  > T_  then  replace 
Ij  Ik  jk  0 


Tjj  in  formulas  (6)  by 


Ij 


TIk  + Tjk* 


(7) 
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If  no  such  k Is  available  Chen  a least  squares  method  can  be  used 
only  with  the  times  available  yielding  formulas  different  from 
equations  (5). 

V.  Conclusions 

A model  for  the  sound  ranging  problem  has  been  formulated  and 
simplified.  A procedure  has  been  presented  for  obtaining  more  ac- 
curate timing  information  from  the  received  microphone  signals.  The 
procedure  described  consisted  of  four  basic  parts: (1)  a pre  filtering 
operation,  (2)  a pair  by  pair  correlation,  (3)  an  elimination  of  in- 
accurate data,  and  (4)  an  overall  least  square  time  fit. 

The  overall  procedure  appears  to  be  feasible  for  operation  but 
not  in  real  time.  Even  not  operating  in  real  time,  the  procedure 
could  be  implemented  at  least  an  order  of  magnitude  faster  than  the 
present  technique.  Also,  operation  in  an  interactive  mode  would  greatly 
minimize  the  possibility  of  presenting  erroneous  timing  information.  The 
prefiltering  operation,  elimination  of  inaccurate  data  and  least  square 
time  fit  parts  can  be  easily  implemented  almost  in  real  time.  The  most 
time  consuming  operation  is  the  pair  by  pair  correlation;  for  example, 
with  six  recorded  signals  we  must  evaluate  fifteen  different  cross  cor- 
relations. These  correlations  would  probably  be  best  done  by  using  the 
product  of  Fourier  Transforms  and  the  inverse  transform.  Existing  soft- 
ware and  hardware  could  be  used  in  a parallel  structure  to  accomplish 
these  cross  correlations  efficiently,  thus  making  the  overall  procedure 
feasible  in  a practical  sense. 
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APPENDIX  A:  Refinement  of  time  differences  estimates 


The  problem  discussed  in  this  appendix  Is  that  of  trying  to  obtain 
the  best  possible  time  difference  estimates  for  a collection  of  time  sig- 
nals (s^(t),  82(c), ...,s(t)}.  Using  s^(t)  as  a reference  we  wish  to  select 
times  t2,t3,...,tN  such  that  the  signals  {s^t) ,s2<t-t2) , . . . ,sN<t-tN) } provide 
the  best  collective  correlation  where  best  will  be  defined  in  terms  of  mini- 
mizing a performance  measure  e. 

The  general  approach  to  this  problem  will  be  to  find  the  T^'s  that 
maximize  the  cross  correlation  function  for  all  pairs  of  signals  s^Ct)  and 
Sj(t)  and  use  these  to  determine  the  t^'s.  The  cross  correlation  function 
Rij  (T)  for  signals  s^t)  and  s^(t)  is  given  by 


VT) 


s1(t)s^ (t-T)dt. 


(A-l) 


Let  X be  the  value  of  T for  which  R^(x)  is  a maximum.  This  specifies  a 

set  of  t ime  values  *^13*  * * * ,*^1N*^23,"^24*  * * * **^2N>  * # * ’ ^N— 1 N ^ 

Now  consider  the  set  of  such  time  values  x|^  for  the  following  set  of  signals 


s|(t)  - s1(t) 

s^(t)  “ si(t-t1)  i = 2, . . . ,N  (A-2) 

The  performance  measure  e that  is  desired  to  be  minimized  is  defined  by 

* " <X12)2  + (T13)2  + •••  + <TiN>2 

+ (T'  )2  + ...  + (T'  )2  (A-3) 


where  the  x*  fs  are  the  values  of  T that  maximize  R±j  (x)  given  by 


r°° 

lij)(T)  " S^(t)8j(t-T)dt. 


(A-4) 


The  TJj's  that  maximize  R^'(T)  given  in  A-4  can  be  found  in  terms  of  the 
t^'s  and  Tij's  by  rewriting  (A-4)  in  terms  of  (A-l).  From  (A-2)  we  can  write 
(A-4)  evaluated  at  as 


Rfj  C ij' 


0° 

Si(t'ti)s(t“tj“Tij)dt 


(A-5) 


By  letting  t - t^  * a in  (A-5)  and  using  (A-l)  we  have 

OO 

Rij1)(TIj)  = si(a)s(a-(T|j+tj-t1))da 


(A-6) 


= Rij (Tij  + W 

Since  maximizes  R^j (T)  we  easily  see  that 

T = t ' + t -t  or  t'  " t + (t  -t  ) (A— 7) 

ij  ij  J i ij  ij  i J 

Substituting  these  values  into  (A-3)  yields  the  following 

e - (t12_C2)  + ^T13~t3)  + •••  + (xiN_tV 

+ (T23-(t2-t3))2  + ...  + (T2N'(t2_tN))2  (A~8) 


+ (XN-l.N+(tN-l"tN)) 


1 <> 


; 


A necessary  and  in  this  case  sufficient  condition,  since  e is  convex, 
is  that 

0 i * 2,3, ... ,N  (A-9) 


9e 

at. 


Taking  the  partial  of  e with  respect  to  the  and  simplifying  gives 


where 


“2  " *’l2  ' T23  ' T24 


T2N] 


“3  ' |T13  + T23  " T34'T35  " " T3N! 


“»  ’ [T1»  + T2N  + Vl,.1 

To  solve  for  the  t^'s  of  (A-10)  we  simply  need  to 

find  the  Inverse  of  the  coefficient  matrix  A(J_  . It  can  be  easily  shown 


that 


aT1  - I 
Vl  N 


2 1 1 ...  1 
1 2 1 1 ...  1 


1 1 


. 1 2 


(A-ll) 


thus  giving  the  t^  vector  that  minimizes  e as 


C2 

V 

t. 

au 

3 

-1 

3 

• 

* Vl 

• 

• 

• 

t,. 

a 

N 

m m 

N 

(A-12) 


“2 

(N-l)  -1  -1  ...  -1 

V 

n ^ 

C2 

a3 

- 

-1  (N-l)  -1  ...  -1 

C3 

= Vi* 

C3 

(A-10) 

> 

-1  -1  ...  -1  (N-l) 

>_ 

>. 

1 
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